function [ results_fnames ] = osl_run_first_level_epoched( oat )

% [ results_fnames ] = osl_run_first_level_epoched( oat )
%
% takes in an OAT, which needs to be setup by calling oat=osl_setup_oat(S), struct
% and runs first level trialwise GLM
%
% This function should normally be called using osl_run_oat(oat);
%
% MWW 2011

global OSLDIR;

first_level=oat.first_level;
source_recon_name='source_recon';

% check if user has asked not to do GLM
if(~isfield(first_level,'doGLM')),
    first_level.doGLM = 1;
end

if(~isfield(first_level,'bc_trialwise')),
    first_level.bc_trialwise=0;
end;

% check if user has asked to jailbreak the trialwise data
if(~isfield(first_level,'save_data'));
    first_level.save_data = 0;
end

if(~isfield(first_level, 'trialwise_directory'));
    first_level.trialwise_directory = [];
end

% set the downsample factor.  If doing a tf, read tf_downsample_factor.  If
% doing time domain, read time_downsample_factor
switch first_level.tf_method
    case 'none' % use time_downsample_factor
        if ~isempty(first_level.time_downsample_factor)
            ds_factor = 1./first_level.time_downsample_factor;
        else
            ds_factor = [];
        end
    otherwise  % use tf_downsample_factor
        if ~isempty(first_level.tf_downsample_factor)
            ds_factor = 1./first_level.tf_downsample_factor;
        else
            ds_factor = [];
        end
end % switch first_level.tf_method

if ~isempty(first_level.time_moving_av_win_size)
    warning('first_level.time_moving_av_win_size is NOT used with epoched data');
end;

contrast_list=first_level.contrast;
for c=1:length(contrast_list),
    contrast_list{c}=contrast_list{c}(:);
end;

clear stats;

if(sum(first_level.bc)>0 && first_level.baseline_timespan(1)>first_level.baseline_timespan(2)),
    error('first_level.baseline_timespan(1)>first_level.baseline_timespan(2)');
end;

if(~isfield(oat.first_level,'design_matrix') && ~isfield(oat.first_level,'design_matrix_summary')),
    error('Design matrix is not specified');
end;

for subi_todo=1:length(first_level.subjects_to_do),   
        
    sub=first_level.subjects_to_do(subi_todo);
    
    %sub/length(first_level.submatfiles),
    
    %% load in the beamformer result for this subject
    source_recon_results_fname = oat.source_recon.results_fnames{sub};
    source_recon_results=osl_load_oat_results(oat,source_recon_results_fname);
    
    % get D and update path in case OAT dir has been moved
    D=source_recon_results.BF.data.D;
    if(~strcmp([oat.source_recon.dirname],D.path)),        
        D=path(D,[oat.source_recon.dirname]);
        D.save;
        D=spm_eeg_load([D.path '/' D.fname]);
        osl_save_oat_results(oat,source_recon_results);    
    end;
    
    %% setup things that are common to all subjects
    if(subi_todo==1),
        
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %% set up mask indices from mask/coords
        if strcmp(source_recon_results.recon_method,'none'), % sensor space analysis
            
            disp('Working in sensor space');
            
            modality='MEG';
            chanindmeg = strmatch(modality, D.chantype);
            first_level_results.mask_indices_in_source_recon=ones(length(chanindmeg),1);
            
        else,
            
            % setup std space brain
            first_level_results.gridstep=source_recon_results.gridstep;
            
            S=[];
            S.lower_level_mask_fname=[oat.source_recon.dirname '/' source_recon_name '_mask'];
            S.current_level_mask_fname=[oat.source_recon.dirname '/' oat.first_level.name '_mask'];
            S.current_level=first_level;
            S.lower_level_mni_coord=source_recon_results.mni_coord;
            S.lower_level_gridstep=source_recon_results.gridstep;
            
            [first_level_results.mask_indices_in_source_recon, first_level_results.mni_coord]=setup_mask_indices(S);
            clear S;
            
        end;
        
        %% set up potential time frequency bases in tf_settings, note that any time or
        % frequency averaging will be done AFTER the time-frequency decomposition
        
        % set the frequency range
        freq_range = first_level.tf_freq_range; % range used for stat testing
        if isempty(freq_range)
            % if not set as a first_level option, default to the range used
            % for source_recon
            freq_range=source_recon_results.source_recon.freq_range;
            
            if isempty(freq_range)
                freq_range=[1 D.fsample/3];
                disp(['WARNING: No frequency range has been specified. Set oat.source_recon.freq_range or oat.first_level.tf_freq_range. Using ' num2str(freq_range) ' Hz.']);                
            end % isempty(freq_range)
        end % if isempty(freq_range)
        
        %source_recon_time_indices = find(D.time>=source_recon_results.source_recon.time_range(1) & D.time<=source_recon_results.source_recon.time_range(2));
        source_recon_time_indices=find(source_recon_results.samples2use);
        source_recon_times = D.time(source_recon_time_indices);
        
        if isempty(first_level.time_range)
            first_level.time_range=source_recon_results.woi;
        end
        
        tf_settings = [];
        tf_settings.tf_method               = first_level.tf_method;
        tf_settings.tf_logtransform         = first_level.tf_logtransform;
        tf_settings.tf_hilbert_freq_res     = first_level.tf_hilbert_freq_res;
        tf_settings.tf_freq_range           = freq_range; % set above
        tf_settings.tf_num_freqs            = first_level.tf_num_freqs;
        tf_settings.tf_multitaper_ncycles   = first_level.tf_multitaper_ncycles;
        tf_settings.raw_times               = source_recon_times;
        tf_settings.ds_factor               = ds_factor;
        tf_settings.tf_morlet_factor        = first_level.tf_morlet_factor;
        tf_settings.tf_multitaper_taper     = first_level.tf_multitaper_taper;
        tf_settings.tf_multitaper_twin      = first_level.tf_multitaper_twin;
        tf_settings.tf_multitaper_freqsmooth= first_level.tf_multitaper_freqsmooth;
        
        out = osl_tf_transform( tf_settings ); % passed without a data matrix, osl_tf_transform just gets container sizes
        
        if isfield(out, 'tf_morlet_basis')
            first_level.tf_morlet_basis = out.tf_morlet_basis;
        end
        
        
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %% setup time and frequency windows/indices for using later
                       
        % setup result time dimension        
        first_level_time_indices  = ( out.tf_times >= first_level.time_range(1) ) & ( out.tf_times <= first_level.time_range(2) );
        first_level_results.times = out.tf_times(first_level_time_indices);                
        if first_level.time_average
            first_level_results.times=mean(first_level_results.times);
            tres=1;
        end
        
        % setup result freq dimension        
        first_level_results.frequencies = out.tf_freqs;
        if first_level.freq_average
            first_level_results.frequencies=mean(first_level_results.frequencies);
        end
        
        % setup result spatial dimension
        Nvoxels_out=length(first_level_results.mask_indices_in_source_recon);        
        if(first_level.space_average)
            Nvoxels_out=1;
            if(isfield(first_level_results,'mni_coord')),
                first_level_results.mni_coord=mean(first_level_results.mni_coord,1);
            end;
        end
        
        clear out
        
        %% create results containers
        
        % get dimensions:
        ntpts                       = length(first_level_results.times);
        nfreqs                      = length(first_level_results.frequencies);
        
        % create cope and stdcope containers:
        first_level_results.stdcope = zeros([Nvoxels_out,ntpts,length(contrast_list),nfreqs]);% num_voxels x  num_timepoints x num_contrasts x num_freqs
        first_level_results.cope    = zeros([Nvoxels_out,ntpts,length(contrast_list),nfreqs]);% num_voxels x num_timepoints x num_contrasts x num_freqs
        first_level_results.pseudo_zstat_var = zeros([Nvoxels_out,1]);
        
    end; % end of if(sub==1)
    
    %% setup the GLM design matrix for this subject
    if(isfield(first_level,'design_matrix')),
        
        first_level.design_matrix;
        
        Ntrials   = size(first_level.design_matrix,1);
        
    else,
        
        % setup trial types
        trialtypes=[];
        triallist=[];

        for i=1:length(source_recon_results.source_recon.conditions), % indexes conditions/triggers within subset
                               
            trigname=source_recon_results.source_recon.conditions{i};
            Ntrialspercond=length(D.pickconditions(trigname)); %% number of trials for this condition
            triallist=[triallist , D.pickconditions(trigname)];           
            trialtypes=[trialtypes ; ones(Ntrialspercond,1).*i];
            
        end;
        
        S3=[];
        S3.trialtypes=trialtypes;
        
        if numel(first_level.design_matrix_summary)==length(oat.source_recon.results_fnames) &&(ischar(first_level.design_matrix_summary{sub})),
            
            S3.Xsummary=first_level.design_matrix_summary{sub};
            
            if(isfield(first_level,'trial_rejects'))
                S3.trial_rejects=first_level.trial_rejects;
            else
                error('Need to set oat.first_level.trial_rejects');
            end;
        else
            S3.Xsummary=first_level.design_matrix_summary;
        end;
        
        first_level.x=setup_beamformer_designmatrix(S3);
        
        Ntrials   = length(triallist);
        
    end;
    
    % precompute some things
    contrast_list=first_level.contrast;
    x=first_level.x;
    pinvxtx=pinv(x'*x);
    pinvx=pinv(x);
    
    %% do everything from now on on one voxel at a time (for the sake of
    %% RAM) in the first level mask space
    
    if first_level.save_data
        % transform-specific parameters
        switch first_level.tf_method
            case 'none'
                trlwise.type = 'time_domain';
            otherwise
                trlwise.type    = first_level.tf_method;
                trlwise.dat     = nan(Ntrials,ntpts,nfreqs); % if not a tf, the last dimenstion will be vestigial
                trlwise.freqs   = first_level_results.frequencies;
                trlwise.times   = first_level_results.times;
        end
        
        % spatial parameters
        trlwise.mask        = read_avw([oat.source_recon.dirname '/' oat.first_level.name '_mask']);
        trlwise.mnicoords   = first_level_results.mni_coord;
        trlwise.gridstep    = first_level_results.gridstep;
    end
    
    %%%%%%%%%%%%%%%
    %% setup channels
    modality='MEG';
    chanindmeg = strmatch(modality, D.chantype);
    chanind = setdiff(chanindmeg, D.badchannels);
    if isempty(chanind)
        error(['No good ' modality ' channels were found.']);
    end
    if strcmp(source_recon_results.recon_method,'none'), % sensor space analysis
        % for sensor space - also include bad sensors sensors 
        first_level_results.mask_indices_in_source_recon=ones(length(chanindmeg),1);
        chanind = chanindmeg;
    end; 
    
    %%%%%%%%%%%%%%%
    %% load in sensor data
    first_level_results.D_sensor_data=D;
    first_level_results.chanind=chanind;
    
    sensor_data=D(chanind, source_recon_time_indices, triallist);
    
    % setup class indices, this allows for multi covs
    NK=numel(source_recon_results.BF.inverse.W.MEG); 
    classchanind=find(strcmp(D.chanlabels,'CLASS'));
    if(isempty(classchanind)),
        error(['No ''CLASS'' chanlabel in: ' D.fname]);
    else
        for kk=1:NK,        
            class_samples_inds{kk} = (D(classchanind,source_recon_time_indices, triallist)==kk);
        end;   
    end;
    
    %%%%%%%%%%%%%%%%
    %% now do space averaging
    if(first_level.space_average),
        
        dat=zeros(Ntrials,length(source_recon_time_indices));
        num=0;
        
        for indind=1:length(first_level_results.mask_indices_in_source_recon), % indexes brain space
            if strcmp(source_recon_results.recon_method,'none'), % sensor space analysis
                weights=zeros(1,length(chanind));
                % for sensor space - just compute for all sensors
                weights(indind)=1;
            else
                weights=source_recon_results.BF.inverse.W.MEG{first_level_results.mask_indices_in_source_recon(indind)};
                first_level_results.pseudo_zstat_var(indind)=weights*weights';
                weights=weights/sqrt(first_level_results.pseudo_zstat_var(indind)); % to compute pseudo z-stats
            end;
            
            if sum(isnan(squash(weights)))==0,
                for tri=1:Ntrials, % indexes trials
                    dat(tri,:)=dat(tri,:)+weights*sensor_data(:, :, tri);
                end;
                num=num+1;
            end;
        end;
        dat=dat/num;
        
    end;            
    
    %% loop over voxels
    hilbert_ds_seed=[];
    ft_progress('init', 'etf');
    
    for indind=1:Nvoxels_out, % indexes brain space
        
        ft_progress(indind/Nvoxels_out);
        
        if(~first_level.space_average),
            dat=nan(Ntrials,length(source_recon_time_indices));

            for kk=1:NK,                                
                if strcmp(source_recon_results.recon_method,'none'), % sensor space analysis
                    weights{kk}=zeros(1,length(chanind));
                    % for sensor space - just compute for all sensors
                    weights{kk}(indind)=1;
                else
                    weights{kk}=source_recon_results.BF.inverse.W.MEG{kk}{first_level_results.mask_indices_in_source_recon(indind)};
                end;
                
                if sum(isnan(squash(weights{kk})))==0,

                    for tri=1:Ntrials, % indexes trials
                        timeinds=find(class_samples_inds{kk}(1,:, tri)); % time indices for class kk
                        if ~isempty(timeinds),
                            dat(tri,timeinds)=weights{kk}*sensor_data(:,timeinds,tri);
                        end;
                    end;
                end;
            end;
        end;
        
        if sum(isnan(squash(dat)))~=0,
            
            first_level_results.cope(indind,:,:,:)=nan;
            first_level_results.stdcope(indind,:,:,:)=nan;
            
        else
            
            %% Transform the data - either downsampling in time domain, or a TF transform
            
            tf_settings = [];
            tf_settings.tf_method               = first_level.tf_method;
            tf_settings.tf_logtransform         = first_level.tf_logtransform;
            tf_settings.tf_hilbert_freq_res     = first_level.tf_hilbert_freq_res;
            tf_settings.tf_freq_range           = freq_range; % set above
            tf_settings.tf_num_freqs            = first_level.tf_num_freqs;
            tf_settings.tf_multitaper_ncycles   = first_level.tf_multitaper_ncycles;
            tf_settings.raw_times               = source_recon_times;
            tf_settings.ds_factor               = ds_factor;
            tf_settings.tf_morlet_factor        = first_level.tf_morlet_factor;
            tf_settings.tf_multitaper_taper     = first_level.tf_multitaper_taper;
            tf_settings.tf_multitaper_twin      = first_level.tf_multitaper_twin;
            tf_settings.tf_multitaper_freqsmooth= first_level.tf_multitaper_freqsmooth;

            if isfield(first_level,'tf_morlet_basis')
                tf_settings.tf_morlet_basis = first_level.tf_morlet_basis;
            end
                        
            out = osl_tf_transform( tf_settings , dat ); % do the transformation
            dat = out.dattf;
            
            if first_level.save_data
                trlwise.dat(indind,:,:,:) = dat;
            end % if first_level.save_data
            
            if(first_level.bc_trialwise),
                baseline_time_indices=first_level_results.times<first_level.baseline_timespan(2) & first_level_results.times>first_level.baseline_timespan(1);
                
                if(sum(baseline_time_indices)>0)
                    % normalise by noise in baseline window:
                    for tri=1:size(dat,1), % indexes trials
                        for f = 1 : size(dat,3),
                            dat(tri,:,f)=(dat(tri,:,f)-mean(dat(tri,baseline_time_indices,f)));
                        end;
                    end;
                else
                    warning('Could not baseline correct as there are no timepoints in the specified baseline period');
                end;
            end;
            
            %% now do windowing over time
            dat=dat(:,first_level_time_indices,:);
            
            %% now do averaging over time and frequency
            if(first_level.time_average)
                dat=mean(dat,2);
            end;
            
            if(first_level.freq_average)
                dat=mean(dat,3);
            end;
            
            %% Do GLM
            if first_level.doGLM
                if(size(dat,1)>1)
                    for f = 1 : size(dat,3),
                        for t= 1 : size(dat,2), % indexes time points within trials
                            
                            if(first_level.use_robust_glm),
                                [b,stats2]=robustfit(x,(dat(:,t,f)),'bisquare',4.685,'off'); % fits GLM robustly, Y=XB+E
                                for c=1:length(contrast_list),
                                    
                                    if(strcmp(oat.first_level.cope_type,'coape')),
                                        first_level_results.cope(indind,t,c,f)=contrast_list{c}'*abs(b);
                                    elseif(strcmp(oat.first_level.cope_type,'cope')),
                                        first_level_results.cope(indind,t,c,f)=contrast_list{c}'*b;
                                    elseif(strcmp(oat.first_level.cope_type,'acope')),
                                        first_level_results.cope(indind,t,c,f)=abs(contrast_list{c}'*b);
                                    else
                                        error('Invalid cope_type');
                                    end;
                                    first_level_results.stdcope(indind,t,c,f)=sqrt(contrast_list{c}'*stats2.covb*contrast_list{c});
                                end;
                                pe(t,:)=b; % num_timepoints x  num_pes
                            else
                                [copeout, varcopeout, coapeout, dof,peout]=glm_fast_for_meg((dat(:,t,f)),x,pinvxtx,pinvx,contrast_list,0);
                                if(strcmp(oat.first_level.cope_type,'coape')),
                                    first_level_results.cope(indind,t,:,f)=coapeout;
                                elseif(strcmp(oat.first_level.cope_type,'cope')),
                                    first_level_results.cope(indind,t,:,f)=copeout;
                                elseif(strcmp(oat.first_level.cope_type,'acope')),
                                    first_level_results.cope(indind,t,:,f)=abs(copeout);
                                else
                                    error('Invalid cope_type');
                                end;
                                
                                first_level_results.stdcope(indind,t,:,f)=sqrt(varcopeout);
                                pe(t,:)=peout; % num_timepoints x  num_pes
                            end;
                            
                        end;
                        
                        %% do bc on the copes
                        if(sum(first_level.bc)>0),
                            baseline_time_indices=first_level_results.times<first_level.baseline_timespan(2) & first_level_results.times>first_level.baseline_timespan(1);
                            
                            for c=1:length(contrast_list),
                                if(first_level.bc(c))
                                    first_level_results.cope(indind,:,c,f)=first_level_results.cope(indind,:,c,f)-mean(first_level_results.cope(indind,baseline_time_indices,c,f),2);
                                end;
                            end;
                            
                        end;
                    end;
                    
                    
                else
                    first_level_results.cope(indind,:,1,:)=dat(1,:,:);
                    first_level_results.stdcope(indind,:,1,:)=1;
                end;
            else
                first_level_results.cope(first_level_results.seed_index,:,:,:)=0;
                first_level_results.stdcope(first_level_results.seed_index,:,:,:)=1e16;
            end;
            
            %if(indind==1)
            %    keyboard;
            %    figure;plot(first_level_results.times,squeeze(first_level_results.cope(1,:,3)));
            %end;
            
        end % if first_level.doGLM
    end % for indind = 1:Nvoxels_out
    
    ft_progress('close');
    
    
    %% save results
    first_level_results.x=x;
    first_level_results.contrasts=first_level.contrast;
    first_level_results.source_recon_results_fname=source_recon_results_fname;
    first_level_results.S=first_level;
    first_level_results.bf_S=source_recon_results.source_recon;
    first_level_results.level=1;
    %first_level_results.source_recon_results_fname=source_recon_results.fname;
    first_level_results.name=first_level.name;
    first_level_results.recon_method=source_recon_results.recon_method;
    first_level_results.subject_name=source_recon_results.subject_name;
    
    % save results to disk
    if first_level.save_data
        if ~isempty(first_level.trialwise_directory)
            trlWise_fname=fullfile(first_level.rawTF_directory, [first_level_results.subject_name '_trlWise']);
        else
            first_level_results.fname=[ first_level_results.subject_name '_' first_level_results.name ];
            trlWise_fname=fullfile(oat.source_recon.dirname, [first_level_results.subject_name '_trlWise']);
        end
        disp(['Saving trialwise data in file ' oat.source_recon.dirname '/' first_level_results.fname]);
        save(trlWise_fname, 'trlwise', '-v7.3');
        clear trlwise
    end % if first_level.save_data
    
    if first_level.doGLM
        first_level_results.fname=[ first_level_results.subject_name '_' first_level_results.name ];
        disp(['Saving statistics in file ' oat.source_recon.dirname '/' first_level_results.fname]);
        osl_save_oat_results( oat, first_level_results );
        results_fnames{sub}=first_level_results.fname;
    end % if first_level.doGLM
    
    if (strcmp(source_recon_results.recon_method,'none') | first_level.space_average),
    else
        disp(['To create niftii files from this use a call to osl_save_nii_stats']);
    end;
end;


